rm(list = ls())

setwd("C:/Users/ym1293/Desktop/Magiya_et_al_Data_and_analysis")

library(readxl)
library(stargazer)
library(survminer)
library(survival)
library(rms)


vali<-read_xlsx("governorships_final_2022_03_01.xlsx", sheet = 1, guess_max = 10000)


#######################################################
#Recoding

vali$Not_right_censored<-(vali$RightCensoredTotal * -1) +1#  Surv takes 0 as right-censored 
vali$vdem_avg<-(vali$vdem_begin+vali$vdem_end)/2#calculates the average of the v-dem liberal democracy score
vali$vdem_elc_avg<-(vali$vdem_elc_start+vali$vdem_elc_end)/2#calculates the average of the v-dem electoral democracy score


median(vali$vdem_avg,na.rm=T)
min(vali$vdem_avg,na.rm=T)
max(vali$vdem_avg,na.rm=T)

median(vali$vdem_elc_avg,na.rm=T)
min(vali$vdem_elc_avg,na.rm=T)
max(vali$vdem_elc_avg,na.rm=T)


#######################################################
#Standardizing v-dem variables

mean.vdem<-mean(vali$vdem_avg,na.rm=T)
sd.vdem<-sd(vali$vdem_avg,na.rm=T)
mean.vdem.elc<-mean(vali$vdem_elc_avg,na.rm=T)
sd.vdem.elc<-sd(vali$vdem_elc_avg,na.rm=T)
vali$vdem_avg<-(vali$vdem_avg-mean.vdem)/sd.vdem#standardizing
vali$vdem_elc_avg<-(vali$vdem_elc_avg-mean.vdem.elc)/sd.vdem.elc#standardizing
#######################################################

#######################################################
#Recoding

vali$gdp_avg<-(vali$GDP_USD_year_end+vali$GDP_USD_year_start)/2#takes average of the local gdp variable
vali$PopPerProvince<-vali$PopPerProvince*100 #converts to percentage
vali$Distance<-vali$Distance/1000 # converts to km


vali$tenure63<-as.numeric(vali$tenure63)
vali$tenure64<-as.numeric(vali$tenure64)
vali$tenure65<-as.numeric(vali$tenure65)
vali$tenure66<-as.numeric(vali$tenure66)
vali$tenure67<-as.numeric(vali$tenure67)
vali$tenure68<-as.numeric(vali$tenure68)
vali$tenure69<-as.numeric(vali$tenure69)
vali$tenure70<-as.numeric(vali$tenure70)
vali$tenure71<-as.numeric(vali$tenure71)
#######################################################


#######################################################
#Subseting for the dataset that has information on the KIM dataset
vali.subset.kim<-subset(vali,YearEnd>1974)
vali.subset.kim<-subset(vali.subset.kim,YearStart<2013)
vali.subset.kim<-subset(vali.subset.kim,total_tenure>0)
vali.subset.kim<-subset(vali.subset.kim,total_tenure>0)
#######################################################


#######################################################
#Recoding

vali.subset.kim$killed.per.month<-vali.subset.kim$Number_Death_KIM/vali.subset.kim$total_tenure#calculates recruit per month
vali.subset.kim$killed.per.month<-log(vali.subset.kim$killed.per.month+1)#log transforms due to skew


erdogan.subset<-subset(vali,tenure63>0|tenure64>0|tenure65>0|tenure66>0|tenure67>0
                       |tenure68>0|tenure69>0|tenure70>0|tenure71>0)#gets subset for Erdogan's period




erdogan.subset$before<-0#creates a variable that equals 1 if the governor was assigned before AKP came into power
erdogan.subset$before[erdogan.subset$tenure62>0]<-1#creates a variable that equals 1 if the governor was assigned before AKP came into power

erdogan.subset$vdem_avg[is.na(erdogan.subset$vdem_avg)]<-erdogan.subset$vdem_begin#calculates average v-dem score for Erdogan's subset

#######################################################
#Subseting for the dataset that has information on the KIM dataset, for Erdogan period's data


erdogan.subset.kim<-subset(erdogan.subset,YearStart>2002)
erdogan.subset.kim<-subset(erdogan.subset.kim,YearStart<2013)
erdogan.subset.kim<-subset(erdogan.subset.kim,total_tenure>0)
#######################################################

#######################################################
#Recoding

erdogan.subset.kim$killed.per.month<-erdogan.subset.kim$Number_Death_KIM/erdogan.subset.kim$total_tenure#calculates number of killed per month
erdogan.subset.kim$killed.per.month<-log(erdogan.subset.kim$killed.per.month+1)#log transforms due to skew


erdogan.subset$Not_right_censored<-(erdogan.subset$RightCensoredAKP * -1) +1#  Surv takes 0 as right-censored 

################################################################
################################################################
################################################################
#####Producing the Models in Table 3############################
################################################################
################################################################
################################################################


#Model 1:
rep.model.dem<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                          vdem_avg
                                        +PopPerProvince
                                        +Distance
                                        +Border
                                        +Sea_Opening
                                        +Slope
                                        +cluster(Province)     
                                        , data = vali)
summary(rep.model.dem)


#Model 2:
rep.model.eth<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                    +province_kurd_50
                                +PopPerProvince
                             +Border
                     +Sea_Opening
                             +Slope
                                +cluster(Province)
                                , data = vali)
summary(rep.model.eth)


#Model 3:
rep.model.killed<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                               +killed.per.month
                             +PopPerProvince
                             +Border
                             +Sea_Opening
                             +Slope
                             +cluster(Province)
                             
                             , data = vali.subset.kim)
summary(rep.model.killed)


#Model 4:
rep.model.elec<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                      +margin_victory
                      +PopPerProvince
                      +Border
                    +Sea_Opening
                    +Slope
                      +cluster(Province)
                      , data = vali)


summary(rep.model.elec)


subset.kim<-vali.subset.kim[complete.cases(vali.subset.kim$killed.per.month),]#subset of complete cases for descriptives
subset.election<-vali[complete.cases(vali$margin_victory),]#subset of complete cases for descriptives



list_means_republic<-list(c("Mean", rep(round(mean(vali$total_tenure, na.rm=T), digits = 2),2)  ,
                            round(mean(subset.kim$total_tenure,na.rm=T), digits=2),
                   round(mean(subset.election$total_tenure, na.rm=T), digits = 2)))


list_sd_republic<-list(c("SD", rep(round(sd(vali$single_tenure, na.rm=T), digits = 2),2),
                         round(sd(subset.kim$total_tenure,na.rm=T), digits=2),
                round(sd(subset.election$single_tenure, na.rm=T), digits = 2)))


#Producing Table 3 with Hazard Ratios
stargazer(rep.model.dem,rep.model.eth,rep.model.killed,rep.model.elec,
          title  = "", apply.coef = exp,  t.auto=F, 
          p.auto=F, report = "vct*",
          covariate.labels = c("Average V-dem Score (1 sd)",
                               "Kurdish Province Dummy",
                                "Kurdish Militant Deaths per Month (logged)",
                                "Govn. Party Elect. Margin in Province",
                               "Province Population pct.",
                                "Distance to Capital (km)",
                                "Border Dummy",
                                  "Sea Opening Dummy",
                                    "Average Slope (degrees)"
          ),
          dep.var.caption  = "Dependent Variable: Governor Tenure Length (Months)",
          add.lines=
            c(list_means_republic,
              list_sd_republic)
)

################################################################
################################################################
################################################################
#####Producing the Models in Table 4############################
################################################################
################################################################
################################################################


erdogan.dem<-coxph(Surv(single_tenure, Not_right_censored) ~ 
                       vdem_avg
                     +PopPerProvince
                      +Distance
                      +Border
                   +Sea_Opening
                   +Slope
                     +before
                     +PopPerPOB
                     +pob_kurd_50
                     +AgeAppoint
                    +gdp_avg
                    +cluster(Province)
                     , data = erdogan.subset)

summary(erdogan.dem)



erdogan.ethnic<-coxph(Surv(single_tenure, Not_right_censored) ~ 
                                 +province_kurd_50
                                 +PopPerProvince
                                  +Border
                                +Sea_Opening
                                +Slope
                                 +before
                                   +PopPerPOB
                                +pob_kurd_50
                                   +AgeAppoint
                                +gdp_avg
                                 +cluster(Province)
                                 , data = erdogan.subset)

summary(erdogan.ethnic)



erdogan.kim<-coxph(Surv(single_tenure, Not_right_censored) ~ 
                           +killed.per.month  
                   +PopPerProvince
                   +Border
                   +Sea_Opening
                   +Slope
                   +before
                   +PopPerPOB
                   +pob_kurd_50
                   +AgeAppoint
                   +gdp_avg
                   +cluster(Province)
                         
                         , data = erdogan.subset.kim)
summary(erdogan.kim)




erdogan.elections<-coxph(Surv(single_tenure, Not_right_censored) ~ 
                               +margin_victory
                         +PopPerProvince
                         +Border
                         +Sea_Opening
                         +Slope
                         +before
                         +gdp_avg
                         +PopPerPOB
                         +pob_kurd_50
                         +AgeAppoint
                         +cluster(Province)
                         , data = erdogan.subset)

summary(erdogan.elections)



subset.election.erdogan<-erdogan.subset[complete.cases(erdogan.subset$margin_victory),]
subset.kim.erdogan<-erdogan.subset.kim[complete.cases(erdogan.subset.kim$margin_victory),]



list_means_erdogan<-list(c("Mean",rep(round(mean(erdogan.subset$single_tenure, na.rm=T), digits = 2),2),
                           round(mean(subset.kim.erdogan$single_tenure,na.rm=T),digits=2),
                   round(mean(subset.election.erdogan$single_tenure, na.rm=T), digits = 2)))

list_sd_erdogan<-list(c("SD", rep(round(sd(erdogan.subset$single_tenure, na.rm=T), digits = 2),2),
                        round(sd(subset.kim.erdogan$single_tenure,na.rm=T),digits=2),
                round(sd(subset.election.erdogan$single_tenure, na.rm=T), digits = 2)))


stargazer(erdogan.dem,erdogan.ethnic,erdogan.kim,erdogan.elections,
          title  = "", apply.coef = exp,  t.auto=F, 
          p.auto=F, report = "vct*",

          covariate.labels = c("Average V-dem Score (1 sd)",
                               "Kurdish Province Dummy",
                               "Kurdish Militant Deaths per Month (logged)",
                              "Govering. Party Electoral Margin in Province",
                              "Province Population pct.",
                              "Distance to capital (km)",
                              "Border Dummy",
                              "Sea Opening Dummy",
                              "Average Slope (degrees)",
                               "Gvrnr Assigned Before AKP Rule",
                              "Province Pop. where Governor Born pct.",
                               "Governor Born in Kurdish Province",
                               "Age Governor Appointed",
                              "Province GDP"
          ),#
          
          dep.var.caption  = "Dependent Variable: Governor Tenure Length (Months)",
          
          add.lines=
            c(list_means_erdogan,
              list_sd_erdogan),
          no.space=TRUE
)         



###################################################
#### For the appendix:
###################################################

################################################################
#Mean Standardizing the V-dem Variables

mean.vdem.begin<-mean(vali$vdem_begin,na.rm=T)
sd.vdem.begin<-sd(vali$vdem_begin,na.rm = T)
vali$vdem_begin_std<-(vali$vdem_begin-mean.vdem.begin)/sd.vdem.begin



mean.vdem.end<-mean(vali$vdem_end,na.rm=T)
sd.vdem.end<-sd(vali$vdem_end,na.rm = T)

vali$vdem_end_std<-(vali$vdem_end-mean.vdem.end)/sd.vdem.end
################################################################


################################################################
################################################################
################################################################
#####Producing the Models in Table A.4##########################
################################################################
################################################################
################################################################



republic.model.vdem.begin<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                 +vdem_begin_std
                               +PopPerProvince
                               +Distance
                               +Border
                               +Sea_Opening
                               +Slope
                               +cluster(Province)
                               , data = vali)
summary(republic.model.vdem.begin)



republic.model.vdem.end<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                              +vdem_end_std
                              +PopPerProvince
                              +Distance
                              +Border
                              +Sea_Opening
                              +Slope
                              
                              +cluster(Province)
                              , data = vali)
summary(republic.model.vdem.end)

stargazer(republic.model.vdem.begin,republic.model.vdem.end,
          title  = "", apply.coef = exp,  t.auto=F, 
          p.auto=F, report = "vct*",
          covariate.labels = c(
            "V-dem Score during the Beginning Year of Assignment",
            "V-dem Score during the Final Year of Assignment",
            "Province Population pct.",
            "Distance to Capital (km)",
            "Border Dummy",
            "Sea Opening Dummy",
            "Average Slope (degrees)"
            
          ),
          
          dep.var.caption  = "Dependent Variable: Governor Tenure Length (Months)",
          
          add.lines=
            c(list_means_republic,
              list_sd_republic),
          no.space=TRUE
          
)

################################################################
################################################################
################################################################
#####Producing the Models in Table A.5##########################
################################################################
################################################################
################################################################


rep.model.vdem.elect<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                       vdem_elc_avg
                                     +PopPerProvince
                                     +Distance
                                     +Border
                                     +Sea_Opening
                                     +Slope
                                     +cluster(Province)     
                                     , data = vali)
summary(rep.model.vdem.elect)



stargazer(rep.model.vdem.elect,
          title  = "", apply.coef = exp,  t.auto=F, 
          p.auto=F, report = "vct*",
          covariate.labels = c(
            "Average V-dem Electoral Democracy Score",
            "Province Population pct.",
            "Distance to Capital (km)",
            "Border Dummy",
            "Sea Opening Dummy",
            "Average Slope (degrees)"
          ),
          
          dep.var.caption  = "Dependent Variable: Governor Tenure Length (Months)",
          no.space=TRUE
)






######################


################################################################
################################################################
################################################################
#####Producing the Models in Table A.6##########################
################################################################
################################################################
################################################################


rep.dem.appndx<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                  vdem_avg
                                +PopPerProvince
                                +Distance
                                +Border
                                +Sea_Opening
                                +Slope
                                +PopPerPOB
                                +pob_kurd_50
                                +AgeAppoint
                                +gdp_avg
                                +cluster(Province)     
                                , data = vali)
summary(rep.dem.appndx)






rep.ethn.appndx<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                               +province_kurd_50
                             +PopPerProvince
                             +Border
                             +Sea_Opening
                             +Slope
                             +PopPerPOB
                             +pob_kurd_50
                             +AgeAppoint
                             +gdp_avg
                             +cluster(Province)
                             , data = vali)
summary(rep.ethn.appndx)


rep.killed.appndx<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                               +killed.per.month
                             +PopPerProvince
                             +Border
                             +Sea_Opening
                             +Slope
                             +PopPerPOB
                             +pob_kurd_50
                             +AgeAppoint
                             +gdp_avg
                             +cluster(Province)
                             , data = vali.subset.kim)
summary(rep.killed.appndx)




rep.elec.appndx<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                  +margin_victory
                                +PopPerProvince
                                +Border
                                +Sea_Opening
                                +Slope
                                +PopPerPOB
                                +pob_kurd_50
                                +AgeAppoint
                                +gdp_avg
                                +cluster(Province)
                                , data = vali)


summary(rep.elec.appndx)



stargazer(rep.dem.appndx,rep.ethn.appndx,
          rep.killed.appndx,rep.elec.appndx,
          #type="html",
          #  out="Hazard Ratios-Republic.html",
          title  = "", apply.coef = exp,  t.auto=F, 
          p.auto=F, report = "vct*",
          #  keep=c("vdem.beg","vdem.end","vdem_avg","province_kurd_50","pob_kurd_50","PopPerProvince",
          #         "military.dummy","Governpartystrentgh"
          #   ),
          covariate.labels = c("Average V-dem Score (1 sd)",
                               #                      "Gvrnr Assigned Before Rule",
                               "Kurdish Province Dummy",
                               "Kurdish Militant Deaths per Month (logged)",
                               "Governing Party Electoral Margin in Province",
                               "Province Population pct.",
                               "Distance to Capital (km)",
                               "Border Dummy",
                               "Sea Opening Dummy",
                               "Average Slope (degrees)",
                               "Province Pop. where Governor Born pct.",
                               "Governor Born in Kurdish Province",
                               "Age Governor Appointed",
                               "Province GDP"
                               
          ),
          dep.var.caption  = "Dependent Variable: Governor Tenure Length (Months)",
          
          add.lines=
            c(list_means_republic,
              list_sd_republic),
          no.space=TRUE
          
)

##################################################################
#Mean Standardizing the V-dem Variables for the AKP/Erdogan Period


mean.vdem.begin.erdogan<-mean(erdogan.subset$vdem_begin,na.rm=T)
sd.vdem.begin.erdogan<-sd(erdogan.subset$vdem_begin,na.rm = T)
erdogan.subset$vdem_begin_std<-(erdogan.subset$vdem_begin-mean.vdem.begin.erdogan)/sd.vdem.begin.erdogan

mean.vdem.end.erdogan<-mean(erdogan.subset$vdem_end,na.rm=T)
sd.vdem.end.erdogan<-sd(erdogan.subset$vdem_end,na.rm = T)
erdogan.subset$vdem_end_std<-(erdogan.subset$vdem_end-mean.vdem.end.erdogan)/sd.vdem.end.erdogan
##################################################################




################################################################
################################################################
################################################################
#####Producing the Models in Table A.7##########################
################################################################
################################################################
################################################################


erdogan.model.vdem.begin<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                   +vdem_begin_std
                                 +PopPerProvince
                                 +Distance
                                 +Border
                                 +Sea_Opening
                                 +Slope
                                +before
                                +PopPerPOB
                                +pob_kurd_50
                                +AgeAppoint
                                +gdp_avg
                                 +cluster(Province)
                                 , data = erdogan.subset)
summary(erdogan.model.vdem.begin)



republic.model.vdem.end<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                 +vdem_end_std
                               +PopPerProvince
                               +Distance
                               +Border
                               +Sea_Opening
                               +Slope
                               +before
                               +PopPerPOB
                               +pob_kurd_50
                               +AgeAppoint
                               +gdp_avg
                               
                               +cluster(Province)
                               , data = erdogan.subset)
summary(republic.model.vdem.end)





stargazer(erdogan.model.vdem.begin,republic.model.vdem.end,
          title  = "", apply.coef = exp,  t.auto=F, 
          p.auto=F, report = "vct*",
          covariate.labels = c(
            "V-dem Score during the Beginning Year of Assignment",
            "V-dem Score during the Final Year of Assignment",
            "Province Population pct.",
            "Distance to Capital (km)",
            "Border Dummy",
            "Sea Opening Dummy",
            "Average Slope (degrees)",
            "Governor Assigned Before Rule",
            "Province Pop. where Governor Born pct.",
            "Governor Born in Kurdish Province",
            "Age Governor Appointed",
            "Province GDP"
          ),
          dep.var.caption  = "Dependent Variable: Governor Tenure Length (Months)",
          add.lines=
            c(list_means_erdogan,
              list_sd_erdogan),
          no.space=TRUE
)

################################################################


################################################################
################################################################
################################################################
#####Producing the Models in Table A.8##########################
################################################################
################################################################
################################################################


erdogan.model.vdem.electoral<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                      vdem_elc_avg
                                    +PopPerProvince
                                    +Distance
                                    +Border
                                    +Sea_Opening
                                    +Slope
                                    +before
                                    +PopPerPOB
                                    +pob_kurd_50
                                    +AgeAppoint
                                    +gdp_avg
                                    +cluster(Province)     
                                    , data = erdogan.subset)
summary(erdogan.model.vdem.electoral)




stargazer(erdogan.model.vdem.electoral,
          #type="html",
          #  out="Hazard Ratios-Republic.html",
          title  = "", apply.coef = exp,  t.auto=F, 
          p.auto=F, report = "vct*",
          covariate.labels = c(
            "Average V-dem Electoral Democracy Score",
            "Province Population pct.",
            "Distance to Capital (km)",
            "Border Dummy",
            "Sea Opening Dummy",
            "Average Slope (degrees)",
            "Governor Assigned Before Rule",
            "Province Pop. where Governor Born pct.",
            "Governor Born in Kurdish Province",
            "Age Governor Appointed",
            "Province GDP"
          ),
          
          dep.var.caption  = "Dependent Variable: Governor Tenure Length (Months)",
          #    add.lines=
          #     c(list_means_erdogan,
          #      list_sd_erdogan),
          no.space=TRUE
)



################################################################
################################################################
################################################################
######Plotting Hazard Ratios-Republic###########################
################################################################
################################################################
################################################################


####Vdem

median.pop.republic<-median(vali$PopPerProvince,na.rm=T)
median.distance.republic<-median(vali$Distance,na.rm=T)
median.border.republic<-median(vali$Border,na.rm=T)
median.sea.opening.republic<-median(vali$Sea_Opening,na.rm = T)
median.slope.republic<-median(vali$Slope,na.rm=T)

minimum.dem.republic<-min(vali$vdem_avg,na.rm=T)
median.dem.republic<-median(vali$vdem_avg,na.rm=T)
maximum.dem.republic<-max(vali$vdem_avg,na.rm=T)



dem_republic_df<- with(vali,
                       data.frame(vdem_avg = c(minimum.dem.republic,median.dem.republic,maximum.dem.republic), 
                                  PopPerProvince = rep(median.pop.republic, 3),
                                  Distance=rep(median.distance.republic,3),
                                  Border=rep(median.border.republic,3),
                                  Sea_Opening=rep(median.sea.opening.republic, 3),
                                  Slope=rep(median.slope.republic,3)
                       )
)

republic.model.democracy.2<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                    vdem_avg
                                  +PopPerProvince
                                  +Distance
                                  +Border
                                  +Sea_Opening
                                  +Slope
                                  +cluster(Province)     
                                  , data = vali)
summary(republic.model.democracy.2)


fit.republic.dem <- survfit(republic.model.democracy.2, newdata = dem_republic_df)




library("ggplot2")
require("tidyr")
ncurve <- ncol(fit.republic.dem$surv)
ssurv <- tidyr::gather(as.data.frame(fit.republic.dem$surv), "strata", "surv")
ssurv$time <- rep(fit.republic.dem$time, ncurve )

vdem_hazard_ratio<-ggplot(data = ssurv) + geom_step(aes(time, surv, color = strata))+
  xlab("Months After Assignment")+ 
  ylab("Probability of Staying in Post")+
  scale_color_manual(labels = c("Min. (0.05)", "Median (0.25)","Max. (0.56)"),
                     values=c("blue","red","green"))+
  labs(color = "V-dem Score")+
  xlim(0, 50)+theme_bw()




####Kurdish Province Dummy

dem_ethnic_df<- with(vali,
                     data.frame(province_kurd_50 = c(0,1), 
                                PopPerProvince = rep(median.pop.republic, 2),
                                Border=rep(median.border.republic,2),
                                Sea_Opening=rep(median.sea.opening.republic, 2),
                                Slope=rep(median.slope.republic,2)
                     )
)

republic.model.ethnic.2<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                 province_kurd_50
                               +PopPerProvince
                               +Border
                               +Sea_Opening
                               +Slope
                               +cluster(Province)     
                               , data = vali)
summary(republic.model.ethnic.2)


fit.ethnic.rep <- survfit(republic.model.ethnic.2, newdata = dem_ethnic_df)


ncurve <- ncol(fit.ethnic.rep$surv)
ssurv <- tidyr::gather(as.data.frame(fit.ethnic.rep$surv), "strata", "surv")
ssurv$time <- rep(fit.ethnic.rep$time, ncurve )

ethnic_hazard_ratio<-ggplot(data = ssurv) + geom_step(aes(time, surv, color = strata))+
  xlab("Months After Assignment")+ 
  ylab("Probability of Staying in Post")+
  scale_color_manual(labels = c("No", "Yes"),
                     values=c("blue","red"))+
  labs(color = "Kurdish Majority\nProvince")+
  xlim(0, 50)+theme_bw()




####KIM Deaths


max.death<-max(vali.subset.kim$killed.per.month,na.rm=T)

dem_death_df<- with(vali,
                    data.frame(killed.per.month = c(0,1,max.death), 
                               PopPerProvince = rep(median.pop.republic, 3),
                               Border=rep(median.border.republic,3),
                               Sea_Opening=rep(median.sea.opening.republic, 3),
                               Slope=rep(median.slope.republic,3)
                    )
)

republic.model.death.2<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                killed.per.month
                              +PopPerProvince
                              +Border
                              +Sea_Opening
                              +Slope
                              +cluster(Province)     
                              , data = vali.subset.kim)
summary(republic.model.death.2)


fit.death.rep <- survfit(republic.model.death.2, newdata = dem_death_df)


ncurve <- ncol(fit.death.rep$surv)
ssurv <- tidyr::gather(as.data.frame(fit.death.rep$surv), "strata", "surv")
ssurv$time <- rep(fit.death.rep$time, ncurve )

death_hazard_ratio<-ggplot(data = ssurv) + geom_step(aes(time, surv, color = strata))+
  xlab("Months After Assignment")+ 
  ylab("Probability of Staying in Post")+
  scale_color_manual(labels = c("0", "1","Max. (3.4)"),
                     values=c("blue","red","green"))+
  labs(color = "Militant\nDeaths\nper Month")+
  xlim(0, 50)+theme_bw()




###Electoral Margin

mean.margin<-mean(vali$margin_victory,na.rm=T)
sd.margin<-sd(vali$margin_victory,na.rm=T)

one.sd.above<-mean.margin+sd.margin
one.sd.below<-mean.margin-sd.margin


margin_df<- with(vali,
                 data.frame(margin_victory = c(one.sd.below,mean.margin,one.sd.above), 
                            PopPerProvince = rep(median.pop.republic, 3),
                            Border=rep(median.border.republic,3),
                            Sea_Opening=rep(median.sea.opening.republic, 3),
                            Slope=rep(median.slope.republic,3)
                 )
)

republic.model.margin.2<-coxph(Surv(total_tenure, Not_right_censored) ~ 
                                 margin_victory
                               +PopPerProvince
                               +Border
                               +Sea_Opening
                               +Slope
                               +cluster(Province)     
                               , data = vali)
summary(republic.model.margin.2)




fit.margin.rep <- survfit(republic.model.margin.2, newdata = margin_df)


ncurve <- ncol(fit.margin.rep$surv)
ssurv <- tidyr::gather(as.data.frame(fit.margin.rep$surv), "strata", "surv")
ssurv$time <- rep(fit.margin.rep$time, ncurve )

margin_hazard_ratio<-ggplot(data = ssurv) + geom_step(aes(time, surv, color = strata))+
  xlab("Months After Assignment")+ 
  ylab("Probability of Staying in Post")+
  scale_color_manual(labels = c("1sd below (-8%)", "Mean (14.3%)","1sd above (%36.7)"),
                     values=c("blue","red","green"))+
  labs(color = "Electoral Margin\nin Province")+
  xlim(0, 50)+theme_bw()




library(patchwork)

republic.graph<-(vdem_hazard_ratio +ethnic_hazard_ratio)/(death_hazard_ratio+margin_hazard_ratio)

#ggsave(republic.graph,file="C:/Users/ym1293/Desktop/hazard_ratio_republic.jpg", height=19.42, width=20.2, units = "cm", dpi=250)




################################################################
################################################################
################################################################
######Plotting Hazard Ratios-Erdogan############################
################################################################
################################################################
################################################################

median.pop.erdogan<-median(erdogan.subset$PopPerProvince,na.rm=T)
median.distance.erdogan<-median(erdogan.subset$Distance,na.rm=T)
median.border.erdogan<-median(erdogan.subset$Border,na.rm=T)
median.sea.opening.erdogan<-median(erdogan.subset$Sea_Opening,na.rm = T)
median.slope.erdogan<-median(erdogan.subset$Slope,na.rm=T)
median.before<-median(erdogan.subset$before,na.rm=T)
median.popperpob<-median(erdogan.subset$PopPerPOB,na.rm=T)
median.pobkurd<-median(erdogan.subset$pob_kurd_50,na.rm=T)
median.age<-median(erdogan.subset$AgeAppoint,na.rm = T)
median.gdp<-median(erdogan.subset$gdp_avg,na.rm=T)



minimum.dem.erdogan<-min(erdogan.subset$vdem_avg,na.rm=T)
median.dem.erdogan<-median(erdogan.subset$vdem_avg,na.rm=T)
maximum.dem.erdogan<-max(erdogan.subset$vdem_avg,na.rm=T)



###Vdem

dem_erdogan_df<- with(erdogan.subset,
                       data.frame(vdem_avg = c(minimum.dem.erdogan,median.dem.erdogan,maximum.dem.erdogan), 
                                  PopPerProvince = rep(median.pop.republic, 3),
                                  Distance=rep(median.distance.republic,3),
                                  Border=rep(median.border.republic,3),
                                  Sea_Opening=rep(median.sea.opening.republic, 3),
                                  Slope=rep(median.slope.republic,3),
                                  before=rep(median.before,3),
                                  PopPerPOB=rep(median.popperpob,3),
                                  pob_kurd_50=rep(median.pobkurd,3),
                                  AgeAppoint=rep(median.age,3),
                                  gdp_avg=rep(median.gdp,3)
                                  
                       )
)



erdogan.model.democracy.2<-coxph(Surv(single_tenure, Not_right_censored) ~ 
                                    vdem_avg
                                  +PopPerProvince
                                  +Distance
                                  +Border
                                  +Sea_Opening
                                  +Slope
                                 +before
                                 +PopPerPOB
                                 +pob_kurd_50
                                 +AgeAppoint
                                 +gdp_avg
                                    +cluster(Province)     
                                  , data = erdogan.subset)
summary(erdogan.model.democracy.2)


fit.erdogan.dem <- survfit(erdogan.model.democracy.2, newdata = dem_erdogan_df)




library("ggplot2")
require("tidyr")
ncurve <- ncol(fit.erdogan.dem$surv)
ssurv <- tidyr::gather(as.data.frame(fit.erdogan.dem$surv), "strata", "surv")
ssurv$time <- rep(fit.erdogan.dem$time, ncurve )


vdem_hazard_ratio_erdogan<-ggplot(data = ssurv) + geom_step(aes(time, surv, color = strata))+
  xlab("Months After Assignment")+ 
  ylab("Probability of Staying in Post")+
  scale_color_manual(labels = c("Min. (0.1)", "Median (0.42)","Max. (0.56)"),
                     values=c("blue","red","green"))+
  labs(color = "V-dem Score")+
  xlim(0, 50)+theme_bw()



###Kurdish Province Dummy


erdogan_ethnic_df<- with(erdogan.subset,
                     data.frame(province_kurd_50 = c(0,1), 
                                PopPerProvince = rep(median.pop.republic, 2),
                                Border=rep(median.border.republic,2),
                                Sea_Opening=rep(median.sea.opening.republic, 2),
                                Slope=rep(median.slope.republic,2),
                                before=rep(median.before,2),
                                PopPerPOB=rep(median.popperpob,2),
                                pob_kurd_50=rep(median.pobkurd,2),
                                AgeAppoint=rep(median.age,2),
                                gdp_avg=rep(median.gdp,2)
                     )
)

erdogan.model.ethnic.2<-coxph(Surv(single_tenure, Not_right_censored) ~ 
                                 province_kurd_50
                               +PopPerProvince
                               +Border
                               +Sea_Opening
                               +Slope
                              +before
                              +PopPerPOB
                              +pob_kurd_50
                              +AgeAppoint
                              +gdp_avg
                                 +cluster(Province)     
                               , data = erdogan.subset)
summary(erdogan.model.ethnic.2)


fit.ethnic.erdogan <- survfit(erdogan.model.ethnic.2, newdata = erdogan_ethnic_df)


ncurve <- ncol(fit.ethnic.erdogan$surv)
ssurv <- tidyr::gather(as.data.frame(fit.ethnic.erdogan$surv), "strata", "surv")
ssurv$time <- rep(fit.ethnic.erdogan$time, ncurve )

ethnic_hazard_ratio_erdogan<-ggplot(data = ssurv) + geom_step(aes(time, surv, color = strata))+
  xlab("Months After Assignment")+ 
  ylab("Probability of Staying in Post")+
  scale_color_manual(labels = c("No", "Yes"),
                     values=c("blue","red"))+
  labs(color = "Kurdish Majority\nProvince")+
  xlim(0, 50)+theme_bw()



#### KIM Deaths


max.death.erdogan<-max(erdogan.subset.kim$killed.per.month,na.rm=T)


erdogan_death_df<- with(erdogan.subset.kim,
                    data.frame(killed.per.month = c(0,1,max.death.erdogan), 
                               PopPerProvince = rep(median.pop.republic, 3),
                               Border=rep(median.border.republic,3),
                               Sea_Opening=rep(median.sea.opening.republic, 3),
                               Slope=rep(median.slope.republic,3),
                               PopPerPOB=rep(median.popperpob,3),
                               pob_kurd_50=rep(median.pobkurd,3),
                               AgeAppoint=rep(median.age,3),
                               gdp_avg=rep(median.gdp,3)
                    )
)

erdogan.model.death.2<-coxph(Surv(single_tenure, Not_right_censored) ~ 
                                killed.per.month
                              +PopPerProvince
                              +Border
                              +Sea_Opening
                              +Slope
                             +PopPerPOB
                             +pob_kurd_50
                             +AgeAppoint
                             +gdp_avg
                                +cluster(Province)     
                              , data = erdogan.subset.kim)
summary(erdogan.model.death.2)


fit.death.rep.erdogan <- survfit(erdogan.model.death.2, newdata = erdogan_death_df)


ncurve <- ncol(fit.death.rep.erdogan$surv)
ssurv <- tidyr::gather(as.data.frame(fit.death.rep.erdogan$surv), "strata", "surv")
ssurv$time <- rep(fit.death.rep.erdogan$time, ncurve )

death_hazard_ratio_erdogan<-ggplot(data = ssurv) + geom_step(aes(time, surv, color = strata))+
  xlab("Months After Assignment")+ 
  ylab("Probability of Staying in Post")+
  scale_color_manual(labels = c("0", "1","Max. (2.1)"),
                     values=c("blue","red","green"))+
  labs(color = "Militant\nDeaths\nper Month")+
  xlim(0, 50)+theme_bw()


ggsave(death_hazard_ratio_erdogan,file="C:/Users/yusuf/My Drive/Ottoman and Republican Bureaucrats/StatisticalAnalysis/yusuf_analysis/New_Plots_Summer_2022/death_hazard_ratio_erdogan.jpg", height=19.42, width=20.2, units = "cm", dpi=250)


###Electoral Margin

mean.margin.erdogan<-mean(erdogan.subset$margin_victory,na.rm=T)
sd.margin.erdogan<-sd(erdogan.subset$margin_victory,na.rm=T)
one.sd.above.erdogan<-mean.margin.erdogan+sd.margin.erdogan
one.sd.below.erdogan<-mean.margin.erdogan-sd.margin.erdogan



margin_df_erdogan<- with(erdogan.subset,
                 data.frame(margin_victory = c(one.sd.below.erdogan,mean.margin.erdogan,one.sd.above.erdogan), 
                            PopPerProvince = rep(median.pop.republic, 3),
                            Border=rep(median.border.republic,3),
                            Sea_Opening=rep(median.sea.opening.republic, 3),
                            Slope=rep(median.slope.republic,3),
                            PopPerPOB=rep(median.popperpob,3),
                            pob_kurd_50=rep(median.pobkurd,3),
                            AgeAppoint=rep(median.age,3),
                            gdp_avg=rep(median.gdp,3)
                 )
)

erdogan.model.margin.2<-coxph(Surv(single_tenure, Not_right_censored) ~ 
                                +margin_victory
                              +PopPerProvince
                              +Border
                              +Sea_Opening
                              +Slope
                              +gdp_avg
                              +PopPerPOB
                              +pob_kurd_50
                              +AgeAppoint
                              +cluster(Province)  
                               , data = erdogan.subset)
summary(erdogan.model.margin.2)




fit.margin.erdogan <- survfit(erdogan.model.margin.2, newdata = margin_df_erdogan)


ncurve <- ncol(fit.margin.erdogan$surv)
ssurv <- tidyr::gather(as.data.frame(fit.margin.erdogan$surv), "strata", "surv")
ssurv$time <- rep(fit.margin.erdogan$time, ncurve )

margin_hazard_ratio_erdogan<-ggplot(data = ssurv) + geom_step(aes(time, surv, color = strata))+
  xlab("Months After Assignment")+ 
  ylab("Probability of Staying in Post")+
  scale_color_manual(labels = c("1sd below (-9.5%)", "Mean (17%)","1sd above (43.5%)"),
                     values=c("blue","red","green"))+
  labs(color = "Electoral Margin\nin Province")+
  xlim(0, 50)+theme_bw()

library(patchwork)

erdogan.graph<-(vdem_hazard_ratio_erdogan +ethnic_hazard_ratio_erdogan)/(death_hazard_ratio_erdogan+margin_hazard_ratio_erdogan)

ggsave(erdogan.graph,file="C:/Users/ym1293/Desktop/hazard_ratio_erdogan.jpg", height=19.42, width=20.2, units = "cm", dpi=250)



